Parte 1: Segmentación de Carreteras en Imágenes Aéreas de Alta Resolución: Desarrollo y Evaluación de un Método Computacional¶
Cargamos las imágenes
import os
from skimage.io import imread
import matplotlib.pyplot as plt
# Directorios
base_path = r"C:\Users\norae\piva_enrique\Materiales1\roads"
sat_path = os.path.join(base_path, "sat")
gt_path = os.path.join(base_path, "gt")
# Listar archivos .tif o .tiff
sat_files = sorted([f for f in os.listdir(sat_path) if f.endswith(".tif") or f.endswith(".tiff")])
gt_files = sorted([f for f in os.listdir(gt_path) if f.endswith(".tif") or f.endswith(".tiff")])
# Comprobar correspondencia entre archivos
assert len(sat_files) == len(gt_files), "Cantidad desigual de imágenes y máscaras"
for i in range(len(sat_files)):
assert sat_files[i].split('.')[0] == gt_files[i].split('.')[0], f"No coinciden: {sat_files[i]} y {gt_files[i]}"
# Inicializar listas
imagenes = []
mascaras = []
# Cargar todas las imágenes y máscaras
for img_file, mask_file in tqdm(zip(sat_files, gt_files), total=len(sat_files), desc="Cargando imágenes"):
img = imread(os.path.join(sat_path, img_file))
mask = imread(os.path.join(gt_path, mask_file))
assert img.shape[:2] == mask.shape[:2], f"Dimensiones diferentes: {img_file}"
imagenes.append(img)
mascaras.append(mask)
print("Número total de imágenes cargadas:", len(imagenes))
print("Tamaño de la primera imagen:", imagenes[0].shape)
print("Tamaño de la primera máscara:", mascaras[0].shape)
# Visualizar ejemplo
plt.figure(figsize=(10, 5))
plt.subplot(1, 2, 1)
plt.imshow(imagenes[0])
plt.title("Imagen satelital")
plt.axis("off")
plt.subplot(1, 2, 2)
plt.imshow(mascaras[0], cmap='gray')
plt.title("Máscara binaria")
plt.axis("off")
plt.show()
Cargando imágenes: 100%|██████████| 20/20 [00:00<00:00, 47.11it/s]
Número total de imágenes cargadas: 20 Tamaño de la primera imagen: (1500, 1500, 3) Tamaño de la primera máscara: (1500, 1500)
Extraemos y visualizamos características
import numpy as np
from skimage.color import rgb2hsv, rgb2lab
from skimage.filters import threshold_local
from scipy.ndimage import gaussian_gradient_magnitude
def extract_features(img):
img = img.astype(np.float32) / 255.0 # Normalización
# Espacios de color
img_hsv = rgb2hsv(img)
img_lab = rgb2lab(img)
# Canales RGB y HSV
r = img[:, :, 0].flatten()
g = img[:, :, 1].flatten()
b = img[:, :, 2].flatten()
h = img_hsv[:, :, 0].flatten()
s = img_hsv[:, :, 1].flatten()
b_lab = img_lab[:, :, 2].flatten()
# Binarizaciones directas
h_bin = (img_hsv[:, :, 0] > 0.2).astype(np.uint8).flatten()
s_bin = (img_hsv[:, :, 1] > 0.3).astype(np.uint8).flatten()
b_bin = (img_lab[:, :, 2] > 20).astype(np.uint8).flatten() # Ajuste porque Lab es distinto en skimage
# Gradientes
mag_h = gaussian_gradient_magnitude(img_hsv[:, :, 0], sigma=1).flatten()
mag_s = gaussian_gradient_magnitude(img_hsv[:, :, 1], sigma=1).flatten()
# Umbralización adaptativa usando skimage
def adaptive_thresh(channel, block_size):
local_thresh = threshold_local(channel, block_size, method='gaussian')
return (channel > local_thresh).astype(np.uint8).flatten()
adapt5 = adaptive_thresh(img[:, :, 0], 5)
adapt15 = adaptive_thresh(img[:, :, 0], 15)
adapt33 = adaptive_thresh(img[:, :, 0], 33)
features = np.stack([
r, g, b, h, s, b_lab,
h_bin, s_bin, b_bin,
mag_h, mag_s,
adapt5, adapt15, adapt33
], axis=1)
return features
import numpy as np
import matplotlib.pyplot as plt
from scipy.ndimage import gaussian_gradient_magnitude
from skimage.color import rgb2hsv, rgb2lab
from skimage.filters import threshold_local
# Imagen normalizada
img = imagenes[0].astype(np.float32) / 255.0
height, width, _ = img.shape
# Espacios de color
img_hsv = rgb2hsv(img)
img_lab = rgb2lab(img)
# Características base
r = img[:, :, 0]
g = img[:, :, 1]
b = img[:, :, 2]
h = img_hsv[:, :, 0]
s = img_hsv[:, :, 1]
b_lab = img_lab[:, :, 2] # Puede tener valores negativos, pero seguimos con él como canal discriminativo
# Gradientes
mag_h = gaussian_gradient_magnitude(h, sigma=1)
mag_s = gaussian_gradient_magnitude(s, sigma=1)
# Umbralizaciones
umbral_h_vals = [0.2, 0.4, 0.6]
umbral_s_vals = [0.2, 0.4, 0.6]
umbral_b_vals = [5, 20, 40] # Ajustados a valores típicos de Lab en skimage
# Adaptativos con skimage
def adaptive_threshold(channel, block_size):
threshold = threshold_local(channel, block_size=block_size, method='gaussian')
return (channel > threshold).astype(np.uint8)
adapt_sizes = [5, 15, 33]
adapt_thresh_r = [adaptive_threshold(r, bs) for bs in adapt_sizes]
# Panel de visualización
features_sets = [
("Canal R", [r]),
("Canal G", [g]),
("Canal B", [b]),
("Canal H", [h]),
("Canal S", [s]),
("Canal *b (Lab)", [b_lab]),
("Umbral H", [(h > t).astype(np.uint8) for t in umbral_h_vals]),
("Umbral S", [(s > t).astype(np.uint8) for t in umbral_s_vals]),
("Umbral *b", [(b_lab > t).astype(np.uint8) for t in umbral_b_vals]),
("Gradiente Mag H", [mag_h]),
("Gradiente Mag S", [mag_s]),
("Adaptativo R", adapt_thresh_r)
]
# Mostrar
cols = 3
rows = len(features_sets)
plt.figure(figsize=(cols * 4, rows * 4))
idx = 1
for name, variants in features_sets:
for i, feat in enumerate(variants):
plt.subplot(rows, cols, idx)
plt.imshow(feat, cmap='gray')
plt.title(f"{name} ({i+1})")
plt.axis('off')
idx += 1
plt.tight_layout()
plt.show()
Características seleccionadas:
Canales R, G, B y H: te dan información de color directa y de tono (H), útil porque la carretera tiene una coloración bastante uniforme.
*Canal b del Lab: separa bien zonas con diferente componente azul-amarillo, y resulta discriminativo en asfalto vs. entorno.
*Umbrales binarios sobre S, b y H: permiten resaltar regiones características (por ejemplo, asfalto tiende a baja saturación y valores intermedios de H).
Gradiente de H: resalta bordes relevantes donde cambia el tono, útil para separar límites de la carretera.
Umbral adaptativo sobre R: adapta el umbral localmente, permitiendo mejorar la segmentación incluso con iluminación desigual.
Características descartadas:
Canal V de HSV: no se incluye porque la luminosidad (V) varía mucho con sombras e iluminación y puede ser sensible al ruido.
Gradiente de S: aporta menos que el de H, pues los cambios de saturación en carretera no siempre delimitan bordes útiles.
Otros umbrales (como H > 0.2 o 0.6): se ve que con 0.4 separa mejor la carretera, así que solo uso ese.
Umbralizaciones múltiples o adaptativas en G o B: tras observar visualmente, no aportaban mejora significativa.
Ya hemos visto que características nos pueden ser útiles, pues ahora usamos solo esas
import numpy as np
from skimage.color import rgb2hsv, rgb2lab
from scipy.ndimage import gaussian_gradient_magnitude
from skimage.filters import threshold_local
def extract_selected_features(img):
# Asegurar formato float en [0, 1]
img = img.astype(np.float32) / 255.0
height, width = img.shape[:2]
# Espacios de color
img_hsv = rgb2hsv(img)
img_lab = rgb2lab(img)
# Características continuas
r = img[:, :, 0].flatten()
g = img[:, :, 1].flatten()
b = img[:, :, 2].flatten()
h = img_hsv[:, :, 0].flatten()
b_lab = img_lab[:, :, 2].flatten()
# Umbralizaciones binarias óptimas
s_bin = (img_hsv[:, :, 1] > 0.2).astype(np.uint8).flatten()
b_bin = (img_lab[:, :, 2] > 20).astype(np.uint8).flatten()
h_bin = (img_hsv[:, :, 0] > 0.4).astype(np.uint8).flatten()
# Gradiente sobre canal H
mag_h = gaussian_gradient_magnitude(img_hsv[:, :, 0], sigma=1).flatten()
# Umbral adaptativo sobre canal R
r_channel = img[:, :, 0]
thresh_adapt15 = threshold_local(r_channel, block_size=15, method='gaussian')
adapt15 = (r_channel > thresh_adapt15).astype(np.uint8).flatten()
# Apilar todas las características
features = np.stack([
r, g, b, h, b_lab,
s_bin, b_bin,
h_bin,
mag_h, adapt15
], axis=1)
return features, height, width
from skimage.io import imread
import os
# Rutas
base_path = r"C:\Users\norae\piva_enrique\Materiales1\roads"
img_path = os.path.join(base_path, "sat", "10078675_15.tiff")
mask_path = os.path.join(base_path, "gt", "10078675_15.tif")
# Cargar imagen y máscara
image = imread(img_path)
mask = imread(mask_path)
# Extraer características y dimensiones
X, h, w = extract_selected_features(image)
y = (mask.flatten() > 128).astype(np.uint8)
# Mostrar dimensiones
print("Tamaño de la imagen:", (h, w))
print("X shape:", X.shape)
print("y shape:", y.shape)
print("Proporción de carretera:", round(np.mean(y), 4))
assert X.shape[0] == y.shape[0], "Las dimensiones de X y y no coinciden"
Tamaño de la imagen: (1500, 1500) X shape: (2250000, 10) y shape: (2250000,) Proporción de carretera: 0.0188
Dividimos en conjunto de entrenamiento y test
import numpy as np
from skimage.transform import resize
# Redimensionar a (1000, 1000)
target_size = (1000, 1000) # (height, width)
pixels_per_image = target_size[0] * target_size[1]
num_total = len(imagenes)
num_test = 2
num_val = 2
num_train = num_total - num_test - num_val
features_selected_list = []
mascaras_resized = []
# Redimensionar imágenes y máscaras + extraer características
for img, mask in zip(imagenes, mascaras):
img_resized = resize(img, target_size, preserve_range=True, anti_aliasing=True).astype(np.uint8)
mask_resized = resize(mask, target_size, order=0, preserve_range=True, anti_aliasing=False).astype(np.uint8)
mascaras_resized.append(mask_resized)
feats, _, _ = extract_selected_features(img_resized)
features_selected_list.append(feats)
# Convertir a arrays
gt_imgs = np.stack(mascaras_resized, axis=0)
y = (gt_imgs.reshape(-1) > 128).astype(np.uint8)
features_selected = np.concatenate(features_selected_list, axis=0)
# División corregida
X_train = features_selected[:pixels_per_image * num_train, :]
X_val = features_selected[pixels_per_image * num_train : pixels_per_image * (num_train + num_val), :]
X_test = features_selected[pixels_per_image * (num_train + num_val):, :]
y_train = y[:pixels_per_image * num_train]
y_val = y[pixels_per_image * num_train : pixels_per_image * (num_train + num_val)]
y_test = y[pixels_per_image * (num_train + num_val):]
# Confirmación
print("features shape:", features_selected.shape)
print("X_train:", X_train.shape)
print("X_val: ", X_val.shape)
print("X_test: ", X_test.shape)
print("y_train:", y_train.shape)
print("y_val: ", y_val.shape)
print("y_test: ", y_test.shape)
features shape: (20000000, 10) X_train: (16000000, 10) X_val: (2000000, 10) X_test: (2000000, 10) y_train: (16000000,) y_val: (2000000,) y_test: (2000000,)
Definición y entrenamiento de modelos
Modelo 1: RNA
from tensorflow import keras
from tensorflow.keras.layers import Dense, Dropout
from tensorflow.keras.optimizers import Adam
from tensorflow.keras.metrics import Accuracy, AUC
def create_model(learning_rate=0.001, hidden_units=64, dropout_rate=0.2):
model = keras.Sequential([
Dense(hidden_units, activation='relu', input_shape=(X_train.shape[1],)),
Dropout(dropout_rate),
Dense(hidden_units, activation='relu'),
Dropout(dropout_rate),
Dense(1, activation='sigmoid') # Capa de salida binaria
])
model.compile(
optimizer=Adam(learning_rate=learning_rate),
loss='binary_crossentropy',
metrics=[
'accuracy',
AUC(name='auc')
]
)
return model
from tensorflow.keras.callbacks import EarlyStopping
from sklearn.utils import class_weight
import numpy as np
# Calcular pesos de clase automáticamente en base a los datos de entrenamiento
class_weights = class_weight.compute_class_weight(
class_weight='balanced',
classes=np.unique(y_train),
y=y_train
)
# Convertir a diccionario como espera Keras
class_weights_dict = dict(enumerate(class_weights))
# Crear el modelo (con varias métricas incluidas si has definido eso en create_model)
model = create_model(learning_rate=0.001)
# Hiperparámetros
batch_size = 1024
epochs = 10
# Entrenamiento con EarlyStopping y class weights
history = model.fit(
X_train, y_train,
validation_data=(X_val, y_val),
batch_size=batch_size,
epochs=epochs,
class_weight=class_weights_dict, # Aplica pesos de clase aquí
callbacks=[
EarlyStopping(
patience=3,
restore_best_weights=True,
monitor='val_loss'
)
],
verbose=1
)
c:\Users\norae\AppData\Local\Programs\Python\Python312\Lib\site-packages\keras\src\layers\core\dense.py:87: UserWarning: Do not pass an `input_shape`/`input_dim` argument to a layer. When using Sequential models, prefer using an `Input(shape)` object as the first layer in the model instead. super().__init__(activity_regularizer=activity_regularizer, **kwargs)
Epoch 1/10 15625/15625 ━━━━━━━━━━━━━━━━━━━━ 63s 4ms/step - accuracy: 0.8559 - auc: 0.9022 - loss: 0.3865 - val_accuracy: 0.8251 - val_auc: 0.8941 - val_loss: 0.4228 Epoch 2/10 15625/15625 ━━━━━━━━━━━━━━━━━━━━ 58s 4ms/step - accuracy: 0.8708 - auc: 0.9161 - loss: 0.3558 - val_accuracy: 0.8256 - val_auc: 0.8973 - val_loss: 0.4155 Epoch 3/10 15625/15625 ━━━━━━━━━━━━━━━━━━━━ 60s 4ms/step - accuracy: 0.8720 - auc: 0.9172 - loss: 0.3540 - val_accuracy: 0.8145 - val_auc: 0.8968 - val_loss: 0.4474 Epoch 4/10 15625/15625 ━━━━━━━━━━━━━━━━━━━━ 61s 4ms/step - accuracy: 0.8726 - auc: 0.9176 - loss: 0.3530 - val_accuracy: 0.8227 - val_auc: 0.8988 - val_loss: 0.4168 Epoch 5/10 15625/15625 ━━━━━━━━━━━━━━━━━━━━ 63s 4ms/step - accuracy: 0.8727 - auc: 0.9183 - loss: 0.3521 - val_accuracy: 0.8230 - val_auc: 0.8995 - val_loss: 0.4344
Hacemos las predicciones con el test
from sklearn.metrics import classification_report, confusion_matrix, roc_curve, roc_auc_score
import matplotlib.pyplot as plt
import numpy as np
from tensorflow.keras.models import load_model
# Predecir probabilidades sobre X_test
y_pred_probs = model.predict(X_test)
y_pred = (y_pred_probs > 0.5).astype(int).flatten()
# Reporte de clasificación
print("Classification Report:")
print(classification_report(y_test, y_pred, digits=4))
# Matriz de confusión
cm = confusion_matrix(y_test, y_pred)
plt.figure(figsize=(6, 5))
plt.imshow(cm, interpolation='nearest', cmap=plt.cm.Blues)
plt.title("Matriz de confusión - Conjunto de Test")
plt.colorbar()
tick_marks = np.arange(2)
clases = ["No carretera", "Carretera"]
plt.xticks(tick_marks, clases)
plt.yticks(tick_marks, clases)
# Anotar cada celda con su valor
thresh = cm.max() / 2.
for i in range(cm.shape[0]):
for j in range(cm.shape[1]):
plt.text(j, i, format(cm[i, j], 'd'),
ha="center", va="center",
color="white" if cm[i, j] > thresh else "black")
plt.ylabel("Etiqueta verdadera")
plt.xlabel("Etiqueta predicha")
plt.tight_layout()
plt.show()
# Curva ROC
fpr, tpr, _ = roc_curve(y_test, y_pred_probs)
auc = roc_auc_score(y_test, y_pred_probs)
plt.figure(figsize=(7, 6))
plt.plot(fpr, tpr, label=f"AUC = {auc:.4f}", color='blue')
plt.plot([0, 1], [0, 1], linestyle="--", color='gray')
plt.title("Curva ROC")
plt.xlabel("Tasa de Falsos Positivos (FPR)")
plt.ylabel("Tasa de Verdaderos Positivos (TPR)")
plt.legend(loc="lower right")
plt.grid(True)
plt.tight_layout()
plt.show()
62500/62500 ━━━━━━━━━━━━━━━━━━━━ 71s 1ms/step Classification Report: precision recall f1-score support 0 0.9934 0.8692 0.9272 1915769 1 0.2260 0.8683 0.3586 84231 accuracy 0.8692 2000000 macro avg 0.6097 0.8688 0.6429 2000000 weighted avg 0.9611 0.8692 0.9032 2000000
Visualizamos las predicciones
import matplotlib.pyplot as plt
# Dimensiones conocidas
h, w = 1000, 1000
# Separar las predicciones binarias por imagen
preds_img1 = y_pred[:h * w].reshape((h, w))
preds_img2 = y_pred[h * w:2 * h * w].reshape((h, w))
# Separar también las máscaras reales
true_mask1 = y_test[:h * w].reshape((h, w))
true_mask2 = y_test[h * w:2 * h * w].reshape((h, w))
# Mostrar resultados
fig, axs = plt.subplots(2, 2, figsize=(12, 12))
axs[0, 0].imshow(true_mask1, cmap='gray')
axs[0, 0].set_title("Máscara real - Imagen 1")
axs[0, 0].axis('off')
axs[0, 1].imshow(preds_img1, cmap='gray')
axs[0, 1].set_title("Predicción - Imagen 1")
axs[0, 1].axis('off')
axs[1, 0].imshow(true_mask2, cmap='gray')
axs[1, 0].set_title("Máscara real - Imagen 2")
axs[1, 0].axis('off')
axs[1, 1].imshow(preds_img2, cmap='gray')
axs[1, 1].set_title("Predicción - Imagen 2")
axs[1, 1].axis('off')
plt.tight_layout()
plt.show()
Hacemos un post procesado
Aplicamos apertura, ya que con erosion primero reducimos el grosor y alrededores de la carretera, eliminando así ruido y después recuperamos tamaño original con dilatación.
import matplotlib.pyplot as plt
import numpy as np
from skimage.morphology import remove_small_objects, erosion, dilation, disk
def postprocesar_apertura2(mask_pred_binaria, min_size=20, erosion_iter=2, selem_radius=1):
"""
Postprocesado morfológico: apertura compuesta por erosión iterada + dilatación,
y eliminación de objetos pequeños.
"""
selem = disk(selem_radius)
mask_bool = (mask_pred_binaria > 0)
# Eliminar objetos pequeños
cleaned = remove_small_objects(mask_bool, min_size=min_size)
# Erosión iterativa
eroded = cleaned
for _ in range(erosion_iter):
eroded = erosion(eroded, selem)
# Dilatación única para recuperar conectividad
dilated = dilation(eroded, selem)
return dilated.astype(np.uint8)
img_shape = (1000, 1000)
num_test_images = 2
pixels_per_image = img_shape[0] * img_shape[1]
for i in range(num_test_images):
start = i * pixels_per_image
end = (i + 1) * pixels_per_image
# Reconstruir máscaras
mask_real = y_test[start:end].reshape(img_shape)
mask_pred = y_pred[start:end].reshape(img_shape)
# Postprocesado
post_mask = postprocesar_apertura2(mask_pred, min_size=20, erosion_iter=2, selem_radius=1)
# Imagen RGB original
img_rgb = imagenes[18 + i]
# Visualización
plt.figure(figsize=(15, 5))
plt.suptitle(f"Visualización postprocesada (Apertura) - Imagen Test #{i+1}", fontsize=16)
plt.subplot(1, 3, 1)
plt.imshow(img_rgb)
plt.title("Imagen original")
plt.axis("off")
plt.subplot(1, 3, 2)
plt.imshow(mask_real, cmap='gray')
plt.title("Máscara real")
plt.axis("off")
plt.subplot(1, 3, 3)
plt.imshow(post_mask, cmap='gray')
plt.title("Predicción postprocesada")
plt.axis("off")
plt.tight_layout()
plt.show()
Modelo 2: Random Forest, con todas las características
from skimage.transform import resize
import numpy as np
target_size = (1000, 1000) # (height, width)
imagenes_resized = []
mascaras_resized = []
for img, mask in zip(imagenes, mascaras):
# Redimensionar imagen con anti-aliasing
resized_img = resize(img, target_size, preserve_range=True, anti_aliasing=True).astype(np.uint8)
# Redimensionar máscara sin suavizado ni interpolación continua
resized_mask = resize(mask, target_size, order=0, preserve_range=True, anti_aliasing=False).astype(np.uint8)
imagenes_resized.append(resized_img)
mascaras_resized.append(resized_mask)
# Reemplazo de las listas originales
imagenes = imagenes_resized
mascaras = mascaras_resized
import numpy as np
features_selected_list = []
y_flat_list = []
# Imágenes 0–17 para entrenamiento (18 imágenes)
for img, mask in zip(imagenes[:18], mascaras[:18]):
feats, _, _ = extract_selected_features(img)
y_flat = (mask.flatten() > 128).astype(np.uint8)
features_selected_list.append(feats)
y_flat_list.append(y_flat)
# Concatenar entrenamiento
X_train_full = np.concatenate(features_selected_list, axis=0)
y_train_full = np.concatenate(y_flat_list, axis=0)
# Test: imágenes 18 y 19
X_test = []
y_test = []
for i in range(18, 20):
feats, _, _ = extract_selected_features(imagenes[i])
y_flat = (mascaras[i].flatten() > 128).astype(np.uint8)
X_test.append(feats)
y_test.append(y_flat)
X_test = np.concatenate(X_test, axis=0)
y_test = np.concatenate(y_test, axis=0)
from sklearn.ensemble import RandomForestClassifier
model_rf = RandomForestClassifier(
n_estimators=100,
max_depth=20,
class_weight='balanced',
random_state=42,
n_jobs=-1,
)
from sklearn.utils import resample
# Extraer igual número de muestras de cada clase
carretera = X_train_full[y_train_full == 1]
no_carretera = X_train_full[y_train_full == 0]
carretera_y = y_train_full[y_train_full == 1]
no_carretera_y = y_train_full[y_train_full == 0]
# Submuestrear negativos (clase mayoritaria)
no_carretera_sample, no_carretera_y_sample = resample(
no_carretera, no_carretera_y,
replace=False,
n_samples=len(carretera), # igualar número
random_state=42
)
# Combinar y barajar
X_bal = np.vstack((carretera, no_carretera_sample))
y_bal = np.hstack((carretera_y, no_carretera_y_sample))
# Entrenar
model_rf.fit(X_bal, y_bal)
RandomForestClassifier(class_weight='balanced', max_depth=20, n_jobs=-1,
random_state=42)In a Jupyter environment, please rerun this cell to show the HTML representation or trust the notebook. On GitHub, the HTML representation is unable to render, please try loading this page with nbviewer.org.
RandomForestClassifier(class_weight='balanced', max_depth=20, n_jobs=-1,
random_state=42)from sklearn.metrics import confusion_matrix, roc_curve, roc_auc_score
import matplotlib.pyplot as plt
import numpy as np
# Predicción de probabilidades y etiquetas binarias
y_proba_rf = model_rf.predict_proba(X_test)[:, 1]
y_pred_rf = (y_proba_rf > 0.5).astype(int)
# Matriz de confusión
cm_rf = confusion_matrix(y_test, y_pred_rf)
plt.figure(figsize=(7, 6))
plt.imshow(cm_rf, interpolation='nearest', cmap=plt.cm.Blues)
plt.title("Matriz de confusión - Conjunto de Test")
plt.colorbar()
clases = ["No carretera", "Carretera"]
tick_marks = np.arange(len(clases))
plt.xticks(tick_marks, clases)
plt.yticks(tick_marks, clases)
# Anotar cada celda
thresh = cm_rf.max() / 2
for i in range(cm_rf.shape[0]):
for j in range(cm_rf.shape[1]):
plt.text(j, i, format(cm_rf[i, j], 'd'),
ha="center", va="center",
color="white" if cm_rf[i, j] > thresh else "black")
plt.ylabel("Etiqueta verdadera")
plt.xlabel("Etiqueta predicha")
plt.tight_layout()
plt.show()
# Curva ROC
fpr, tpr, _ = roc_curve(y_test, y_proba_rf)
auc_rf = roc_auc_score(y_test, y_proba_rf)
plt.figure(figsize=(7, 6))
plt.plot(fpr, tpr, color='blue', label=f'AUC = {auc_rf:.4f}')
plt.plot([0, 1], [0, 1], linestyle='--', color='gray')
plt.title("Curva ROC - Random Forest")
plt.xlabel("Tasa de Falsos Positivos (FPR)")
plt.ylabel("Tasa de Verdaderos Positivos (TPR)")
plt.legend(loc='lower right')
plt.grid(True)
plt.tight_layout()
plt.show()
import matplotlib.pyplot as plt
import numpy as np
# Parámetros
img_shape = (1000, 1000)
num_test_images = 2
pixels_per_image = img_shape[0] * img_shape[1]
for i in range(num_test_images):
start = i * pixels_per_image
end = (i + 1) * pixels_per_image
# Máscara real y predicción binaria para esta imagen
mask_real = y_test[start:end].reshape(img_shape)
pred_bin = y_pred_rf[start:end].reshape(img_shape)
# Usar imágenes 18 y 19
img_rgb = imagenes[18 + i]
# Visualización
plt.figure(figsize=(15, 5))
plt.suptitle(f"Predicciones Random Forest - Imagen Test #{i+1}", fontsize=16)
plt.subplot(1, 3, 1)
plt.imshow(img_rgb)
plt.title("Imagen original")
plt.axis("off")
plt.subplot(1, 3, 2)
plt.imshow(mask_real, cmap='gray')
plt.title("Máscara real")
plt.axis("off")
plt.subplot(1, 3, 3)
plt.imshow(pred_bin, cmap='gray')
plt.title("Predicción binaria RF")
plt.axis("off")
plt.tight_layout()
plt.show()
En este postprocesado solo aplicamos erosión, para ver cómo cambia sin la necesidad de aplicar una dilatación posterior. Saldrán las lineas de las carreteras más finas.
import numpy as np
from skimage.morphology import remove_small_objects, erosion, disk
def postprocesar_erosion(pred_binaria, radio_kernel=1, min_size=50):
"""
Postprocesado morfológico con erosión usando skimage.
Parámetros:
- pred_binaria: array 2D binario (0 o 1)
- radio_kernel: radio del elemento estructurante (disk)
- min_size: tamaño mínimo de objetos que se mantienen
Retorna:
- Máscara binaria erosionada
"""
# Binarizar
mask = (pred_binaria > 0).astype(bool)
# Eliminar objetos pequeños
mask = remove_small_objects(mask, min_size=min_size)
# Crear elemento estructurante
selem = disk(radio_kernel)
# Aplicar erosión
mask_eroded = erosion(mask, selem)
return mask_eroded.astype(np.uint8)
img_shape = (1000, 1000)
num_test_images = 2
pixels_per_image = img_shape[0] * img_shape[1]
for i in range(num_test_images):
start = i * pixels_per_image
end = (i + 1) * pixels_per_image
# Reconstruir máscara binaria
pred_binaria = y_pred_rf[start:end].reshape(img_shape)
# Aplicar solo erosión
post_mask2 = postprocesar_erosion(pred_binaria, radio_kernel=1, min_size=50)
# Visualización
plt.figure(figsize=(10, 5))
plt.suptitle(f"Postprocesado con erosión - Imagen Test #{i+1}", fontsize=14)
plt.subplot(1, 2, 1)
plt.imshow(pred_binaria, cmap='gray')
plt.title("Predicción binaria RF")
plt.axis("off")
plt.subplot(1, 2, 2)
plt.imshow(post_mask2, cmap='gray')
plt.title("Solo erosión (skimage)")
plt.axis("off")
plt.tight_layout()
plt.show()
Ahora quiero probar cogiendo como características solo los canales del modelo RGB
Modelo 3: Random Forest con RGB
import numpy as np
# Entrenamiento: imágenes 0–17 (18 imágenes)
X_train_rgb, y_train_rgb = [], []
for i in range(18): # imágenes 0 a 17
rgb_feats = imagenes[i].reshape(-1, 3)
y_flat = (mascaras[i].flatten() > 128).astype(np.uint8)
X_train_rgb.append(rgb_feats)
y_train_rgb.append(y_flat)
X_train_rgb = np.concatenate(X_train_rgb, axis=0)
y_train_rgb = np.concatenate(y_train_rgb, axis=0)
# Test: imágenes 18 y 19 (2 imágenes)
X_test_rgb, y_test_rgb = [], []
for i in range(18, 20):
rgb_feats = imagenes[i].reshape(-1, 3)
y_flat = (mascaras[i].flatten() > 128).astype(np.uint8)
X_test_rgb.append(rgb_feats)
y_test_rgb.append(y_flat)
X_test_rgb = np.concatenate(X_test_rgb, axis=0)
y_test_rgb = np.concatenate(y_test_rgb, axis=0)
from sklearn.utils import resample
import numpy as np
# Separar clases
X_carretera_rgb = X_train_rgb[y_train_rgb == 1]
y_carretera_rgb = y_train_rgb[y_train_rgb == 1]
X_no_carretera_rgb = X_train_rgb[y_train_rgb == 0]
y_no_carretera_rgb = y_train_rgb[y_train_rgb == 0]
# Submuestrear negativos para igualar número
X_no_carretera_bal_rgb, y_no_carretera_bal_rgb = resample(
X_no_carretera_rgb,
y_no_carretera_rgb,
replace=False,
n_samples=len(y_carretera_rgb),
random_state=42
)
# Combinar clases balanceadas
X_train_bal_rgb = np.vstack((X_carretera_rgb, X_no_carretera_bal_rgb))
y_train_bal_rgb = np.hstack((y_carretera_rgb, y_no_carretera_bal_rgb))
# Mezclar aleatoriamente
indices_rgb = np.random.permutation(len(y_train_bal_rgb))
X_train_bal_rgb = X_train_bal_rgb[indices_rgb]
y_train_bal_rgb = y_train_bal_rgb[indices_rgb]
from sklearn.ensemble import RandomForestClassifier
model_rf_rgb = RandomForestClassifier(
n_estimators=100,
max_depth=20,
class_weight='balanced',
random_state=42,
n_jobs=-1
)
model_rf_rgb.fit(X_train_bal_rgb, y_train_bal_rgb)
RandomForestClassifier(class_weight='balanced', max_depth=20, n_jobs=-1,
random_state=42)In a Jupyter environment, please rerun this cell to show the HTML representation or trust the notebook. On GitHub, the HTML representation is unable to render, please try loading this page with nbviewer.org.
RandomForestClassifier(class_weight='balanced', max_depth=20, n_jobs=-1,
random_state=42)from sklearn.metrics import classification_report, confusion_matrix, roc_curve, roc_auc_score
import matplotlib.pyplot as plt
import numpy as np
# Obtener las probabilidades de predicción
y_proba_rgb = model_rf_rgb.predict_proba(X_test_rgb)[:, 1]
# Convertir a predicción binaria
y_pred_rgb = (y_proba_rgb > 0.5).astype(int)
# Reporte de clasificación
print("Classification Report - RGB Features:")
print(classification_report(y_test_rgb, y_pred_rgb, digits=4))
# Matriz de confusión
cm_rgb = confusion_matrix(y_test_rgb, y_pred_rgb)
plt.figure(figsize=(6, 5))
plt.imshow(cm_rgb, interpolation='nearest', cmap=plt.cm.Blues)
plt.title("Matriz de Confusión - Test RGB")
plt.colorbar()
tick_marks = np.arange(2)
clases = ["No carretera", "Carretera"]
plt.xticks(tick_marks, clases)
plt.yticks(tick_marks, clases)
# Añadir los valores en las celdas
thresh = cm_rgb.max() / 2
for i in range(cm_rgb.shape[0]):
for j in range(cm_rgb.shape[1]):
plt.text(j, i, format(cm_rgb[i, j], 'd'),
ha="center", va="center",
color="white" if cm_rgb[i, j] > thresh else "black")
plt.ylabel("Etiqueta verdadera")
plt.xlabel("Etiqueta predicha")
plt.tight_layout()
plt.show()
# Curva ROC
fpr, tpr, _ = roc_curve(y_test_rgb, y_proba_rgb)
auc = roc_auc_score(y_test_rgb, y_proba_rgb)
plt.figure(figsize=(7, 6))
plt.plot(fpr, tpr, label=f"AUC = {auc:.4f}", color='darkorange')
plt.plot([0, 1], [0, 1], linestyle="--", color='gray')
plt.title("Curva ROC - Random Forest (RGB)")
plt.xlabel("Tasa de Falsos Positivos (FPR)")
plt.ylabel("Tasa de Verdaderos Positivos (TPR)")
plt.legend(loc="lower right")
plt.grid(True)
plt.tight_layout()
plt.show()
Classification Report - RGB Features:
precision recall f1-score support
0 0.9938 0.8607 0.9224 1915769
1 0.2168 0.8772 0.3477 84231
accuracy 0.8614 2000000
macro avg 0.6053 0.8689 0.6351 2000000
weighted avg 0.9610 0.8614 0.8982 2000000
import matplotlib.pyplot as plt
# Parámetros de la imagen
img_shape = (1000, 1000)
num_test_images = 2
pixels_per_image = img_shape[0] * img_shape[1]
for i in range(num_test_images):
start = i * pixels_per_image
end = (i + 1) * pixels_per_image
# Reconstruir máscaras
pred_img = y_pred_rgb[start:end].reshape(img_shape)
mask_real = y_test_rgb[start:end].reshape(img_shape)
img_rgb = imagenes[18 + i]
# Mostrar visualización
plt.figure(figsize=(15, 5))
plt.suptitle(f"Predicciones RF - Características RGB - Imagen Test #{i+1}", fontsize=16)
plt.subplot(1, 3, 1)
plt.imshow(img_rgb)
plt.title("Imagen original")
plt.axis("off")
plt.subplot(1, 3, 2)
plt.imshow(mask_real, cmap='gray')
plt.title("Máscara real")
plt.axis("off")
plt.subplot(1, 3, 3)
plt.imshow(pred_img, cmap='gray')
plt.title("Predicción binaria RF")
plt.axis("off")
plt.tight_layout()
plt.show()
Aquí volvemos a aplicar apertura, ya que antes vimos que da mejores resultados que sólo erosión.
from skimage.morphology import remove_small_objects, opening, disk
import numpy as np
import matplotlib.pyplot as plt
def postprocesar_apertura2(pred_binaria, min_size=300, radio_kernel=1):
"""
Aplica limpieza morfológica: eliminación de objetos pequeños y apertura morfológica.
"""
mask = (pred_binaria > 0).astype(bool)
mask = remove_small_objects(mask, min_size=min_size)
selem = disk(radio_kernel)
mask = opening(mask, selem)
return mask.astype(np.uint8)
# Parámetros
img_shape = (1000, 1000)
num_test_images = 2
pixels_per_image = img_shape[0] * img_shape[1]
for i in range(num_test_images):
start = i * pixels_per_image
end = (i + 1) * pixels_per_image
# Predicción binaria
pred_bin = y_pred_rgb[start:end].reshape(img_shape)
mask_real = y_test_rgb[start:end].reshape(img_shape)
img_rgb = imagenes[18 + i] # imágenes de test
# Postprocesamiento
post_mask3 = postprocesar_apertura2(pred_bin, min_size=300, radio_kernel=1)
# Visualización
plt.figure(figsize=(15, 5))
plt.suptitle(f"Postprocesado RF RGB - Imagen Test #{i+1}", fontsize=16)
plt.subplot(1, 3, 1)
plt.imshow(img_rgb)
plt.title("Imagen original")
plt.axis("off")
plt.subplot(1, 3, 2)
plt.imshow(mask_real, cmap='gray')
plt.title("Máscara real")
plt.axis("off")
plt.subplot(1, 3, 3)
plt.imshow(post_mask3, cmap='gray')
plt.title("Predicción postprocesada")
plt.axis("off")
plt.tight_layout()
plt.show()
Superposición de los resultados sobre las imágenes satelitales
def mostrar_resultados_modelo(img, mask, index_ignorado, title):
img = img.astype(np.uint8)
mask = mask.reshape(img.shape[0], img.shape[1])
overlay = img.copy()
overlay[mask == 1] = [255, 255, 0]
plt.figure(figsize=(12, 4))
plt.subplot(1, 3, 1)
plt.imshow(img)
plt.title("Imagen original")
plt.axis("off")
plt.subplot(1, 3, 2)
plt.imshow(mask, cmap='gray')
plt.title("Máscara predicha")
plt.axis("off")
plt.subplot(1, 3, 3)
plt.imshow(overlay)
plt.title(title)
plt.axis("off")
plt.tight_layout()
plt.show()
mostrar_resultados_modelo(imagenes[19], post_mask1, 0, "Modelo 1")
mostrar_resultados_modelo(imagenes[19], post_mask2, 0, "Modelo 2")
mostrar_resultados_modelo(imagenes[19], post_mask3, 0, "Modelo 3")